function [x1,f1] = VecFun_GoldenSearch(f,param,a,b,tol)

alpha1 = (3-sqrt(5))/2;
alpha2 = (sqrt(5)-1)/2;
d  = b-a;
x1 = a+alpha1*d;
x2 = a+alpha2*d;
f1 = f(x1,param);
f2 = f(x2,param);

d = alpha1*alpha2*d;
while max(d)>tol
  d = d*alpha2;
  % x2 is new upper bound
  I_x2  =   f2<f1;
  if any(I_x2)
      x2(I_x2)  =  x1(I_x2);
      x1(I_x2)  =   x2(I_x2)-d(I_x2);
      f2(I_x2)  =   f1(I_x2);
      f1(I_x2)  =   f(x1(I_x2),param(I_x2,:));
  end
  % x1 is new lower bound
  I_x1      =   ~I_x2;
  if any(I_x1)
      x1(I_x1)  =   x2(I_x1);
      x2(I_x1)  =   x1(I_x1)+d(I_x1);
      f1(I_x1)  =   f2(I_x1);
      f2(I_x1)  =   f(x2(I_x1),param(I_x1,:));
  end
end

% Return the larger of the two
I_max       =   f2>f1;
if any(I_max)
    x1(I_max)   =   x2(I_max);
    f1(I_max)   =   f2(I_max);
end
